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Abstract. We compute spectra of symmetric random matrices defined on graphs 
exhibiting a modular structure. Modules are initially introduced as fully connected sub- 
units of a graph. By contrast, inter-module connectivity is taken to be incomplete. Two 
different types of inter-module connectivity are considered, one where the number of 
intermodule connections per-node diverges, and one where this number remains finite in 
the infinite module-size limit. In the first case, results can be understood as a perturbation 
of a superposition of semicircular spectral densities one would obtain for uncoupled 
modules. In the second case, matters can be more involved, and depend in detail on 
inter- module connectivities. For suitable parameters we even find near-triangular shaped 
spectral densities, similar to those observed in certain scale-free networks, in a system 
of consisting of just two coupled modules. Analytic results are presented for the infinite 
module-size limit; they are well corroborated by numerical simulations. 



1. Introduction 

Appreciation has steadily grown in recent years that theories of networked systems 
provide useful paradigms to understand complex processes in various branches of science 
and technology, including the evolution of the internet and the world-wide web, flows 
of information, power, or traffic, credit, market, or operational risk, food-webs in 
ecosystems, gene-regulation, protein-protein interactions underlying metabolic processes 
or cell signalling, immune system response, information processing in neural networks, 
opinion formation, the adoption of new technologies in societies, or the spread of diseases 
or epidemics, and more (see, e.g. [H [2] for recent reviews). 

Random matrix theory has long been known to constitute a powerful tool to study 
topological properties of the graphs underlying networked systems 0, HI El [lj E]. E.g., 
moments of the spectral density of an adjacency matrix describing a graph, gives complete 
information about the number of walks returning to the originating vertex after a given 
number of steps. 
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Any form of heterogeneity existing in a graph, be it due to a scale-free distribution of 
connectivities [TJ, due to a small- world structure [8], or due to the existence of a modular 
structure in terms of identifiable sub-graphs (9j, or a certain density of specific motifs 
existing in a graph would, therefore, manifest itself in the spectral density of states of 
the corresponding connectivity matrix [3j HI El [lj El [9] . A triangular density of states, for 
instance, has been hailed as a signature of a scale-free connectivity distribution emerging 
from network growth by preferential attachment [5]. 

Spectra of complex networks with heterogeneity arising from scale-free or small-world 
connectivity patterns were obtained through numerical diagonalization [3 [10]; analytical 
results for certain networks with power-law degree distributions exist, based either 
on the effective medium approximation(EMA) [TT], or restricted to the limit of large 
average connectivity [121 13] • A complete solution for general degree distributions with 
finite average connectivity (scale-free or other) has been obtained only recently, both 
for ensembles in the thermodynamic limit [TJ], and for large single instances [15J, the 
latter generalised to non-hermitean matrices in [16]. A highly efficient approximation for 
Poissonian random graphs at moderately large mean connectivity has appeared in [TT] . 
We are not aware of studies other than numerical for the case of networks with small 
world [5 J or explicitly modular structure [9]. 

In the present paper, we introduce a block-structured random matrix model to study 
spectral properties of coupled systems, or of a system with modules. The model can 
be thought of as a network composed of sub-networks, with weighted links within and 
between subnetworks. To simplify matters, we will initially consider the case where 
modules are defined as fully connected subnetworks. However, we note at the outset 
that this case is not substantially different from the one where intra-module connectivity 
is incomplete but with coordination that diverges in the infinite module-size limit. Inter- 
module connectivity is assumed to be incomplete, and two substantially different cases 
will be considered, one where the number of inter-module connections per node diverges, 
and one where this number remains finite in the infinite module-size limit. 

We shall find that results for the first case can be rationalized in terms of perturbed 
superpositions of semi-circular spectral densities characteristic of the individual modules 
when considered on their own, whereas in the second case matters can be more involved, 
and depend in greater detail on inter-module connectivities. For suitable parameters we 
even find near-triangular shaped spectral densities, similar to those observed in certain 
scale-free networks, already for of two coupled modules. 

The remainder of the paper is organized as follows. In Sec. 2 we present an elementary 
variant of a modular system a system consisting of two coupled modules. This section 
mainly serves to introduce the formal definitions, and the method we will be using 
for computing the density of states, viz. the the replica method as proposed in the 
random matrix context by Edwards and Jones [18]. Details of the replica calculation 
for the extensively cross-connected case are given in Sec. 2.1. The fundamentally different 
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case of (two) finitely cross-connected modules is dealt with in Sec. 2.2. In Sec. 3 we 
present a generalisation of the self-consistency equations derived earlier for the general 
M-modules case. Results for the extensively and finitely cross-connected two-modules case 
are presented and discussed in Sees 4.1 and 4.2, respectively. The paper concludes with a 
summary and outlook in Sec. 5. 



2. A Modular System 

We begin by considering the most elementary modular system described by a 2N x 2N 
symmetric block matrix of the form 

in which the symmetric N x N sub- matrices M^ 2 \ and the N x N coupling-sub- 

matrix V are random. We assume that elements of M^ 1 ' and are independently and 
identically Gaussian distributed: M\f ~ 7V(0, Jj/iV), 1 < % < j < N, /i = 1, 2, and that 
V i:j = djKij with Kij rsj JV(o, Jp/c), 1 < i,j < N, while 

p( c h) = ( x - ^) <^,o + fy^i ■ (2) 

Two different limits will be considered: (i) the so-called extensively cross-connected limit 
c — > oo, —> oo, with the ratio c/N either remaining finite or approaching zero in the 
N — > oo-limitjj] (ii) the so-called finitely cross-connected limit, where c is kept constant, 
as iV — > oo. In this case the distribution of inter-module connectivities is Poissonian in 
the thermodynamic limit, with average coordination c. 

We are interested in the spectral density of H, 

2N 

M^) = ^E^-^) > (3) 



fc=l 



more precisely in its average Pat (A) over the random matrix ensemble introduced above. 
Here, the are the eigenvalues of H. We shall use p(A) to denote the (average) spectral 
density in the thermodaynamic limit, 



p(A) = lim PN (X) . (4) 

N— >oo 

The spectral density is computed from the resolvent via 

PnW = RS 2A^ Im ^ lndet [K1 ~ HVl/2 ' (5) 



| The second alternative could more appropriately be referred to as sub-extensive, a distinction we are 
not going to make here for simplicity, as it does not affect the nature of results 
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in which X £ = A — is, and the inverse square root of the determinant is obtained as a 
Gaussian integral. Using and to denote N component vectors, and u = (W 1 -*, u^) 
to denote their concatenation, we get 



Ptv(A) = lim Im-^- /in 



(2ti/i) n 6XP 



[A e I-if]u 



(6) 



where angled brackets on the r.h.s denote an average over connectivities {cy} and weights 
{M^ ] } and {Kij} of the non-vanishing matrix elements. 



The average of the logarithm is evaluated using replica. 
p N (X) = lim 



— Im^- lim - \n(ZZ) 



(7) 



with 



N 



Here a = 1, . . . , n enumerates the replica. 



0=1 



S.i. T/ie Extensively Cross- Connected Case 



The average (Zjv) is easily performed [18], [19]. We have 



(expj^E^-^^ 
^expj.E^"^}) 



exp UE«i 2K ^ (2) «i 2) 



Up to subdominant corrections from diagonal matrix elements this gives 



(z n N ) 



^ 2 



E(^ 1} • «i lh2 



n 



i+ 



J; 



N 



exp 



2c 



«a "ifi "ja 



(9) 



Expanding the exponential in the product (for large c) and re-exponentiating one obtains 



(z n N ) 



n 



(27T/C 



^ exp l-i £ A.(«W • «?> + ni 2 ) ■ W i 2 )) - ^ • -i 1} ) 2 



E(-i 2) • -i 2) ) 2 - § E(-i 1} • ■ *4 2) ) 

a, ft a, b 



(10) 
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N 



E 



„.(m)„m) 

u ia a ib 



H=l,2 



as order parameters and by enforcing their definition via (^-functions. This leads to 

iridic) 



with 



(z n N ) 

G 2 = 



II t ^P{N[Gi + G 2 + G 3 ]} 



^i,a,b 



2tt/N 



J 



(2)n2 



Jl 



06 



ab 



ah 



ab 



G, = In 



1 Z-^ylab Hab + <?a6 9, 

dnl^dn?) 



aft 



/"tT dU a 'dWa ' Z r 



fab) U a " >U \> ^ 



-^U^a, b -2^)U^ 
ab 

- ~lndet(A e I-2?«) - ~ In det(A £ I - 2g( 2 )) 
2 2 



:i2i 



(13) 
(14) 



(15) 



2.1.1. Replica Symmetry — Fixed Point Equations In the large N limit the density of 
states is dominated by the saddle-point contribution to Eq. (fl2l) . Adopting the by now 
well-established assumptions of replica-symmetry, and invariance under rotation in the 
space of replica at the relevant saddle-point [18] 



„(a0 _ _0)i 

Qab - Qd d ab 



dp) _ ;Wr 
- Id °ab 



(16) 



one has 



G — G\ + G*2 4" G3 

' 2 " #(«f) 2 



T 2 

J P (1) (2) 

-fTd <ld 



■ (1) -~(2) (2) 



ln(A £ -2g«) + ln(A £ -2gf; 



(17) 



where terms of order n 2 are omitted. 

Stationarity of G requires that the RS order parameters solve the following fixed point 
equations 



■ 41) 



1 



-Jh { d x) + ~JpQd 



1 



2J2) 



-(2) 



. T 2 a {2) 4- - T A a K 

, J 2Hd + J pQ.d 



2„W 
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Qd 



1 



i\ £ - 2iq 



(i) 



i(2) 



i\ £ - 2iq\ 

These can be combined by eliminating the conjugate variables to give, 

1 (2) 1 



tf 3 



i\ e + JfqV + .P p q d 



(2) ' 



i\ £ + J*q? + Jlq d l 



;is) 



2.1.2. Spectral Density for the c — > oo Case The average spectral density is obtained by 
differentiation w.r.t A via Eq. ([7j). Only terms with explicit A-dependence in G of Eq. (TTJ 
contribute at the saddle point. This gives 



,(A) = -Re 



(19) 



For the limiting cases J p = 0, describing two uncoupled systems, and J\ = J2 = J, 
describing a coupling of two systems with identical statistics of intra-system couplings we 
obtain the following explicit analytical results. 

In the uncoupled case J p = the solution of Eqs. (ITS]) is 



± 



4 7 2 - A 2 



giving 



P(A) 



1 

27 



e(4j 2 - a 2 



4 Jl - A 2 + 



9(4 J 2 - A 2 
2J 2 



4J| - A 2 



(20) 



(21) 



9 P 

i.e. a superposition of two independent Wigner semi-circular densities as expected. 

In the case of two coupled systems which are statistically identical, with J\ = J 2 = J 
the solution of Eqs. f|T8|) is q d = q d = qd with 



Qd 



i\ f 



± 



1 



leading to 



P(A) 



2( J 2 + J 2 ) 2(J 2 + J 2 ) 
6(4(J 2 + J 2 ) - A 2 ^ 



4 (J2 + J2) _ A 



f22l 



4(J 2 + J 2 ) - A 2 



2tt( J 2 + J 2 ) 

i.e. a Wigner semi-circular density, with the radius of the semi-circle given by r 
2 J J 2 + Jl 



(23) 



In the asymmetric case J\ 7^ J2 we solve the fixed point equations numerically. A non-zero 
coupling leads to a smoothing of the superposition of the two independent semi-circles as 
anticipated; see Sec. 0] below. 
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2.2. Two Finitely Cross-Connected Modules 

We now consider a system with the same basic setup, except that we assume the average 
number of cross-connections c to be finite (also in the thermodynamic limit). In this 
case, the analysis is considerably more involved. It combines techniques used in (18] for 
connected and of [19] for sparse random matrices, and more specifically the reformulation 
[13] of the sparse case that allows to proceed to explicit results. In this case the calculations 
can be carried out without restricting the distribution of the non-zero matrix elements of 
the inter-module connections to be Gaussian, and we will in developing the theory not 
make that restriction. 

For the average of the replicated partition function we get 

dwi 1} dwi 2) _ J l ^ , (1) (2) (2) , _ J[_ (1) (Ih2 



w> - n Toi/7; ex p e ^k i} • «?> + • < } ) - ek ; 



7 2 

2 



X>i 2) • -i 2) ) 2 + ^E ( (exp {^E^})k " ') } (24) 



AN , 

a,b 

in analogy to Eq. (191), where ( ... )# represents an average over the It^ distribution, 
which is as yet left open. As in Eq. ([9]), subdominant contributions coming from diagonal 
matrix elements are omitted. 

Decoupling of sites is achieved by introducing order parameters 

tf^^E^^ , M=l,2, (25) 

i 

as in the extensively cross-connected case before, but in addition also the replicated 
densities 



p<W(u) v 



^EII<^-4 M) ) , ^ = 1,2, (26) 



as well as the corresponding conjugate quantities. This allows to express Eq. (I2"lj) as 
an integral over the order parameters and their conjugates, combined with a functional 
integral over the replicated densities and their conjugates, 



(z n N ) = J mvp^w^} J n dq ; b J q / ^{n^+g.+g,]} , (27) 

M fi,a,b ' 



with 



Gi= -|e(^) 2 -|e(^) 2 



G 2 = -iE(?M + ?2 ) ?£ ) )- i /dn(p«(n)p«(n)+p( 2 )( W )p( 2 )( W ) 
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+1 »/n-«.«pU<«: 



06 



Here we have introduced abbreviations of the form dp^\u) = du for integrals 

over densities where appropriate to simplify notation. 



2.2.1. Replica Symmetry & Self- Consistency Equations Eq. (J2"T1) is evaluated by the 
saddle point method. The saddle point for this problem is expected to be replica- 
symmetric and symmetric under rotation in the space of replica as in the extensively 
cross-connected case. In the present context this translates to an ansatz of the form 



(lab — Id °ab , 



%b 



(2* 



for the Edwards Anderson type order parameters [18J and 

exp[-| M 2] 



p^\u) 



ip^\u) 



dvr^Hn- 



Z(io) 



Z(u) 



(29) 



i.e. an uncountably infinite superposition of complex Gaussians (with Reu; > and 
Rea) > 0) for the replicated densities and their conjugates [Hj. Here we have introduced 
the shorthand 



du exp 



-U 



^2tt/uj . 



(30) 



The in the expressions for pV*) are to be determined such that the densities fc^ are 
normalised. 

This ansatz translates path-integrals over the replicated densities p and p into path- 
integrals over the densities 7r and tt, giving 



with now 



[Z n N ) 



Gi ~ n 



2n/N 



exp{N[G 1 + G 2 + G 3 ]} (31) 



4 



,/ ' — rw<i ) 



2 ^2 /„(-•) ^ 



g. 



?? 



r.^(i) (i) . ..(2) (2) 



2 

-£ 



M +nc {fl) /d7r M (^)d7r ( ^(cj)ln 



Z(u)Z(u) 



(32) 



(33) 
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+ B £ P4W (t)/{d#M} t h,^» 
fc=o J lU=i^{v e ) 

Here {d7r^} fc = flLi dTr^(^), while = Eti^, and 
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(34) 



PcM(k) = — — exp[-c w ] 



fc! 



(35) 



is a Poissonian distribution with average (k) = and we have introduced the partition 
functions 



exp 



i\ e -2i$ ) + {u} k )u* 

1/2 



JX e - 2iq { f } + {uj} k 



— - ( w« 2 + a/-u 2 — 2/K tic 



Z 2 (uj,uj', K) = J dudv exp 
The stationarity conditions for tt^(lu) and ti^(uj) then read 



2tt 



vW + fT 2 ' 



(36) 



£37) 



In 



+ a;) 
Z(u)Z(u) 



and 



| dTT^) 



In 



+ a;) 
Z(u)Z(u) 



=°/^w{*t$$} k + » (38) 



with Hi and /z 2 Lagrange multipliers to take the normalisation of and into account. 
The stationarity conditions for the 7r^(u)) are 



; y d7r M (cj)ln 



Z{Cj + c<j) 



EfcP^W /{ d ^Vil° ^j^^ + ^ (40) 
fc >i 7 il = i 



where the cr M are once more Lagrange multipliers to take the normalisation of the (u>) 
into account. 

The stationarity conditions for the qf^ are 



*9d 



/* 2 

2 9d 



(41) 



the corresponding ones for the conjugate variables give 

7 — n ^ I Jfc 



k=0 

oo 



k=o J i 



i\ £ + J^f + {u} k 



(42) 
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in which (. . is defined as an average w.r.t. the Gaussian weight in terms of which 



tt (1) (£) = / d7i i2 \uj')(5{u-Q{u',K))) K (43) 



Z&*> ({£>}*) i s defined. We have used Eq. gT]) to express the in terms of the qf* . 

Following [201 EI], the stationarity conditions for and 7r( 2 )(ci;) are rewritten in a 

form that suggests solving them via a population based algorithm. In the present case we 
get 



and 



in which 



while 



ft^(u) = / d^ 1 \u')(5(u-Q(u',K))) K , (44) 



~ K 2 

n = — , (45) 

Co"' 



E / {d*<"W (u> - «W) (46) 

fc>l c J 



with 



fc-i 



£=1 



In Eqs. ( )43l) . ()44l) . and (|4"6"I) . we have already invested = c to enforce that the it^ 
are normalized. 

For the spectral density we obtain the same formal result as for the extensively cross- 
connected system before, 



1 



p(A) = -Re [g« + gf j . (48) 

The solution of these coupled sets of equations is considerably more involved than for the 
extensively cross-connected systems considered earlier, as it involves solving equations for 
the coupled macroscopic order parameters q^' and q^\ which are themselves expressed 
in terms of averages over self-consistent solutions of a pair of non-linear integral equations 
parameterised by these order parameters. However, we have found that a population 
dynamics in which values of q^ and q^ are iteratively updated using Eq. (T4"2|) by sampling 
from the corresponding populations rapidly converges to the correct solution. 
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3. The Multi-Modular Case 



Finally, we consider systems with M modules, each of size N, mutually cross-connected 
with finite connectivity. Inside modules we may or may not have all-to-all Gaussian 
couplings of variances Jt/N. In addition there are finitely many 0(1) module-to- module 
couplings for each vertex, with average connectivities cf v for couplings between nodes in 
modules \i and v (possibly including the case fj, — v) . Evaluating this case requires a fairly 
straightforward generalisation of the setup developed earlier. Here we only produce the 
fixed point equations and the final expression for the average spectral density. We get 



7T 



=E^T / d^\u){5(u-n(u,K))) (49) 

C J /* 



and 



7>) 



u ) = E ~M k ) I i^h-iS (w - n&) , (so) 

k>l C J 



as well as 



with 



and 



$ = EpAQ {d^h , n W ~- k n (51) 



K 2 

Q(u,K) = — , (52) 



n& = iA. + Jj# ) + (53) 

t=\ 

In Eq. flUD , (. . .)fa, denotes an average over the distribution of couplings connecting 
modules fi and u, and we have the normalisation J2v c M1/ = c M . 

For the spectral density we obtain the (obvious) generalisation of those obtained for the 
two modules case before, viz. 

i M 

P(A) = ^ReE^ } - (54) 
[i — i 

Generalising these to the situation of varying module-sizes = r^N, with r M > 0, is 
straightforward. 
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Figure 1. Spectral densities for J\ — 1, J 2 = 0.1, and J p — 0.5 (left), and J p — 0.75 
(right). 

4. Results 

4-1. Extensively Cross- Connected Systems 

The results we obtain for extensively cross-connected systems can be understood in terms 
of superpositions of Wigner semi-circles one would expect for the spectral densities of 
the modules if they were uncoupled, but smoothed (and broadened) by the interaction. 
We have checked that our results agree perfectly with simulations, but have not included 
results of simulations in the figures for the extensively cross-connected systems below. 

In Fig. 1 we explore the effect of the inter-module coupling strength J p on the spectral 
density of a system of two coupled modules with intra-module coupling strengths J M 
differing by one order of magnitude. The reader is invited to compare the results with 
those expected for in a situation where the two modules were non-interacting, namely a 
simple superposition of semi-circular densities of radii 1.1 \ and 2J 2 , respectively. 

In Fig. 2 we keep the strength of the inter-module couplings but vary the ratio of the two 
intra-module couplings. 

Results for uncoupled modules or for coupled modules with identical intra-module 
coupling statistics are not shown; they were found to be in perfect agreement with the 
simple analytical results presented in Sec. 2.1.2. 

4-2. Finitely Cross-Connected Systems 

The modifications generated by cross-connections that remain finite in number and 
strength for each node in each of the blocks are more pronounced than those created 
by extensive (infinitesimal) cross-connections. We mention just two fairly drastic 
modifications. First, the spectral density of a system of finitely cross-connected modules 
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Figure 2. Spectral densities for Ji = 1, J 2 = 0.25 (left), and J 2 = 0.5 (right), and 
J p = 0.5 in both cases 

does no longer have sharp edges with square-root singularities of the spectral density at 
the edges as in the case of extensively cross-connected systems, where it derives from a 
perturbed superposition of semi-circular densities, but rather tails with decay laws that 
depend on the nature of the distribution of cross-connections. In the present case of 
Poisson distributions of cross-connections we find these tails to exhibit exponential decay. 
Second, even with identical intra-module coupling statistics the spectral density is no 
longer given by a simple semi-circular law as it is for the extensively cross-connected case, 
but rather exhibits marked deviations from the semi-circular law, with details depending 
both on the number and the strengths of the cross-connections, as shown in Fig. 3. 

As the numerics in the present finitely cross-connected case is considerably more 
involved we present results of the analytic theory together with checks against numerical 
diagonalization. Figs 3-4 demonstrate that the analytic results are in excellent agreement 
with numerical simulations, virtually indistinguishable for the parameters and the 
statistics used. 

Fig 3 exhibits spectral densities of a two-module system with identical intra-module 
coupling statistics, and Gaussian cross-connections with a Poissonian degree statistics, of 
average cross-coordinations 2 and 5 respectively. In the first case the spectral density has 
a shape close to triangular, though more rounded at the tip than what is known from 
spectral density of certain scale-free systems. 

Fig. 4 looks at a system of two modules having intra-module connections of different 
strengths, and average cross-coordinations 2; results resemble the corresponding case 
of extensive cross-connectivity, apart from the exponential tails, which are an exclusive 
feature of the finitely cross-connected case. 
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Figure 3. Spectral densities for J\ = J 2 = 1 and J p = \j\fc for c = 2 (left) and 
c = 5 (right). Results of numerical diagonalizations of 500 matrices containing two 
coupled blocks, each of dimension 1000, are shown for comparison (dashed lines); they 
are virtually indistinguishable from results of the analytic theory. 




Figure 4. Spectral density for J\ = 1, J 2 = 0.5 and J p — 1/y/c for c = 2 (full line). 
Results of numerical diagonalizations of 500 matrices containing two coupled blocks, each 
of dimension 1000, are shown for comparison (dashed line), and are once more virtually 
indistinguishable from results of the analytic theory. 



5. Summary and Conclusions 

We have evaluated spectral densities of symmetric matrices describing modular systems. 
Modularity is regarded as one of several routes to create heterogeneity in interacting 
systems. In some biological systems in fact, modularity of interactions appears to be a 
natural consequence of compartmentalization; systems with cellular structure, or sub- 
structures within cells come to mind, where heterogeneity of interaction patterns due 
to modularity of the system would seem to enjoy a greater degree of plausibility than 
heterogeneity as observed in certain scale free systems. 

In any case, whenever large systems with different levels of organization are considered, 
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modularity appears to be a feature to be reckoned with. 

While the results for extensively cross-connected systems can be understood in terms of 
perturbed superpositions of semi-circular spectral densities characteristic of the individual 
modules if in isolation, the case of finitely cross-connected is more sensitive to details of 
the statistics of cross-connections. 

We have evaluated examples only for systems containing two coupled modules, but 
presented the general theory for multi- modular systems in Sec 3. Our derivations and 
results were restricted to the case where entries of the connectivity matrix were sampled 
independently (leading to Poissonian degree distributions of inter-module connections 
in the finitely cross-connected case). Unlike in the case of a single module [22], the case 
where degree distributions of inter-module connections differ from Poissonian distributions 
requires more substantial modifications to the structure of the theory We will deal with 
these, as well as with spectra of matrices corresponding to small-world systems [8] in a 
separate publication [23] . 
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